LES of Temporally Evolving Mixing Layers by an 
Eighth- Order Filter Scheme 

A. Hadjadj*’ a , H.C. Yee b , B. Sjogreen c 

a CORIA UMR 6614 & INS A de Rouen, 76800 St-Etienne du Rouvray, France 
b NASA Ames Research Center, Moffett Field, CA, 94035, USA 
c Lawrence Livermore National Laboratories, Livermore, CA, 94551, USA 


Abstract 

An eighth-order filter method for a wide range of compressible flow speeds (H.C. Yee and 
B. Sjogreen, Proceedings of ICOSAHOM09, June 22-26, 2009, Trondheim, Norway) are 
employed for large eddy simulations (LES) of temporally evolving mixing layers (TML) 
for different convective Mach numbers ( M c ) and Reynolds numbers. The high order fil- 
ter method is designed for accurate and efficient simulations of shock-free compressible 
turbulence, turbulence with shocklets and turbulence with strong shocks with minimum 
tuning of scheme parameters. The value of M c considered is for the TML range from the 
quasi-incompressible regime to the highly compressible supersonic regime. The three main 
characteristics of compressible TML (the self similarity property, compressibility effects 
and the presence of large-scale structure with shocklets for high M c ) are considered for the 
LES study. The LES results using the same scheme parameters for all studied cases agree 
well with experimental results of Barone et al. (2006), and published direct numerical 
simulations (DNS) work of Rogers & Moser (1994) and Pantano & Sarkar (2002). 
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1. Motivation and Objective 

Lor the last decade there has been an increase in the use of computational fluid dynam- 
ics (CLD) in engineering science, not only for fundamental understanding of complex com- 
pressible turbulent physics, but also for the development and design of industrial devices. 
Owing to the recent progress in petascale computing, in tandem with advances in algorithm 
development for accurate direct numerical simulations (DNS) and large eddy simulations 
(LES) of shock free compressible turbulence and turbulence with strong shocks, this type of 
DNS and LES computation has gradually been able to tackle more complex flows physics. 
Advances in flow visualization tools have paved the way to extracting valuable informa- 
tion from the computed results containing hundreds of terabytes of data. Examples include 
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flows through internal propulsive nozzles with shock-wave propagation or sound emission 
from supersonic jets, and mixing and shock/boundary layer interactions. 

In compressible turbulent combustion/nonequilibrium flows, the constructions of nu- 
merical schemes for (a) stable and accurate simulation of turbulence with strong shocks, and 
(b) obtaining correct propagation speed of discontinuities for stiff reacting terms on “coarse 
grids” share one important ingredient - minimization of numerical dissipation while main- 
taining numerical stability. Here “coarse grids” means standard mesh density requirement 
for accurate simulation of typical non-reacting flows. This dual requirement to achieve 
both numerical stability and accuracy with zero or minimal use of numerical dissipation is 
most often conflicting for existing schemes that were designed for non-reacting flows. In 
addition to the minimization of numerical dissipation while maintaining numerical stability 
in compressible turbulence with strong shocks, Yee & Sjogreen, Yee and Yee & Sweby 
[65, 66, 62, 61, 68, 69] discussed a general framework for the design of such schemes. 
Yee & Sjogreen [70], Sjogreen & Yee [53] and Wei et al. [60] and references cited therein 
present their recent progress on the subject. In [73], a short overview of this recent progress 
is given. The discussion addresses three separate yet interwoven types of numerical chal- 
lenges for high speed turbulent reacting flows containing discontinuities. This paper is 
confined to the study of turbulent mixing for non-reacting flows. The study for turbulent 
mixing for reacting flows is planned. 

1.1. Recent Progress in Numerical Methods for Turbulence with Strong Shocks 

The current trends in the containment of numerical dissipation in DNS and LES of tur- 
bulence with shocks are summarized in Yee & Sjogreen and Yee et al. [70, 69, 72]. See the 
cited references for details on these current trends. Before presenting the temporally evolv- 
ing mixing layers (TML) studies, the key ingredients and the performance of the high order 
nonlinear filter schemes with pre-processing and post-processing steps in conjunction with 
the use of a high order non-dissipative spatial base scheme [70, 72] are briefly illustrated 
for two test cases. 

1.1.1. High Order Nonlinear Filter Schemes [51, 68, 70, 72] 

Before the application of a high order non-dissipative spatial base scheme, the pre- 
processing step to improve stability had split inviscid flux derivatives of the governing 
equation(s) in the following three ways, depending on the flow types and the desire for 
rigorous mathematical analysis or physical argument. 

• Entropy splitting of Olsson & Oliger [38] and Yee et al. [64, 65]: The resulting 
form is non-conservative and the derivation is based on entropy norm stability with 
numerical boundary closure for the initial value boundary problem. 

• The system form of the Ducros et al. splitting [12]: This is a conservative splitting 
and the derivation is based on physical arguments. 
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• Tadmor entropy conservation formulation for systems (Sjogreen & Yee [54]): The 
derivation is based on mathematical analysis. It is a generalization of Tadmor’s en- 
tropy formulation to systems and has not been fully tested on complex flows. 

After the application of a non-dissipative high order spatial base scheme on the split 
form of the governing equation(s), to further improve nonlinear stability from the non- 
dissipative spatial base scheme, the post-processing step of Yee & Sjogreen [68, 70], 
Sjogreen & Yee [51] nonlinearly filtered the solution by a dissipative portion of a high order 
shock-capturing scheme with a local flow sensor. These flow sensors provide locations and 
amounts of built-in shock-capturing dissipation that can be further reduced or eliminated. 
For all the computations shown, the Ducros et al. splitting is employed since a conserva- 
tive splitting is more appropriate if one does not know if the subject flow is shock-free or 
turbulence with shocks. Some attributes of the high order filter approach are: 

• Spatial Base Scheme: High order and conservative with high order freestream preser- 
vation metric evaluation for curvilinear grids, (no flux limiter or Riemann solver) 

• Physical Viscosity: Automatically taken into consideration by the base scheme. The 
same order of central differencing for the viscous derivative as the convective flux 
derivatives are used. 

• Efficiency: One Riemann solve per dimension per time step, independent of time dis- 
cretizations (less CPU time and fewer grid points than their standard shock-capturing 
scheme counterparts) 

• Accuracy: Containment of numerical dissipation via local wavelet flow sensor 

• Well-balanced scheme: These nonlinear filter schemes are well-balanced schemes for 
certain chemical reacting flows and problem containing geometric source terms [59] 

• Parallel Algorithm: Suitable for most current supercomputer architectures 

1.2. Sample test Cases Illustrating the Efficiency and Accuracy of High Order Filter 
Schemes 

These filter schemes are efficient, and the total computational cost for a given error tol- 
erance is lower than for standard shock-capturing schemes of the same order. This is of 
importance, for example, in DNS and in flow control optimization to improve aerodynamic 
properties, where the flow simulation must be carried out many times during the optimiza- 
tion loop. The efficiency and accuracy of the schemes for a wide variety of flow problems 
can be found in aforementioned cited references. Here two test cases are illustrated. 

2D Shoek/Vorticity Interaction: Figure 1 shows a comparison of a second-order TVD, 
seventh-order WENO (WEN07), hybrid scheme (switch between eighth-order spatial cen- 
tral scheme and WEN07 using wavelet flow sensor as the switch indicator) and the filter 
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Figure 1: CPU comparison of four shock-capturing schemes. 


scheme WEN07fi (an eighth-order spatial central base scheme and the dissipative por- 
tion of WEN07 and the same wavelet flow sensor to guide where the WEN07 dissipation 
should be applied at the post-possessing nonlinear filter step). A second-order Runge-Kutta 
method was used for the TVD scheme and the classical fourth-order Runge Kutta method 
was used for the rest of the spatial scheme. For this particular simple 2D shock-vorticity 
interaction test case with a simple weak planar shock without structure, WEN07, hybrid, 
and WEN07fi give the same accuracy. However, there is large gain in CPU time by the 
filter scheme for this turbulence-free test case. For turbulence with shocks, there is a more 
beneficial gain both in accuracy and CPU time of the filter schemes over the their standard 
WENO counterparts. 

ID Shock/Turbulence Interaction Problem: This 1-D compressible inviscid ideal gas 
problem is one of the most computed test cases in the literature to assess the capability of a 
shock-capturing scheme in the presence of shock/turbulence interactions. The flow consists 
of a shock at Mach 3 propagating into a sinusoidal density field with initial data given 
by (p L , ul , Pl) = (3.857143, 2.629369, 10.33333) to the left of a shock located at 
x = —4, and (p R , u R , p R ) = (1 + 0.2 sin(5x), 0, 1) to the right of the shock, where p is 
the density, u is the velocity and p is the pressure. The computational domain is [—5, 5] and 
the computation stops at time equal to 1.8. Figure 2 shows the comparison among WEN03, 
WEN05 and WEN07, and their corresponding filter schemes WEN03fi, WEN05fi and 
WEN07fi using a very coarse uniform grid of 200 points with the reference solution. The 
reference solution is obtained with WEN05 using 16000 grid points. WEN05fi required 
at the most 50% of the CPU time of WEN05 if third or fourth-order Runge-Kutta time 
discretization were used. In order for WEN05 to obtain a similar accuracy as WEN05fi, at 
least two times the number of grid points is needed. Moreover, the accuracy of WEN05fi 
is similar to WEN09 (computation not shown). 
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Figure 2: 


1 . 3 . Objective and Outline 

In this paper we report recent studies of LES computations of compressible turbulent 
TML flows using numerical schemes developed by [51, 68, 70] in conjunction with the pre- 
processing step discussed previously. The current research is motivated by the overarching 
goal of developing numerical tools for reliable predictive capability of complex turbulent 
reacting flows, especially for problems including compressibility, heat transfer and real gas 
effects interacting with instabilities, shocks and turbulence. The comparative study among 
WEN07A, WEN05 and WEN07 is reported in [71] with grid refinement studies was the 
first step in determining the suitable order of filter schemes to be used for the current physics 
based study. The LES filtering issue in the presence of shocks [47] is not addressed. The 
paper is organized as follows: Numerical methods are given in Section 2. The subgrid 
models for the compressible Navier-Stokes equations are given briefly in Section 3 and 
Appendix A. Results are then presented and discussed in Section 4. 

2. Numerical Method 

This section summarizes the numerical methods to be used for the turbulent TML study. 
The numerical methods solve the split form of the inviscid flux derivatives according to the 
pre-processing step. The discussion is broken up into two subsections. 
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2.1. Original High Order Filter Method 


For turbulence with shocks, instead of solely relying on very high order high-resolution 
shock-capturing methods for accuracy, the filter schemes [63, 64, 51, 67, 68] take advantage 
of the effectiveness of the nonlinear dissipation contained in good shock-capturing schemes 
as stabilizing mechanisms at locations where needed. Such a filter method consists of two 
steps: a full time step using a spatially high-order non-dissipative base scheme, followed 
by a post-processing filter step. The post-processing filter step consists of the products 
of wavelet-based flow sensors and nonlinear numerical dissipations. The flow sensor is 
used in an adaptive procedure to analyze the computed flow data and indicate the location 
and type of built-in numerical dissipation that can be eliminated or further reduced. The 
nonlinear dissipative portion of a high-resolution shock-capturing scheme can be any TVD, 
MUSCL, ENO, or WENO scheme. By design, the flow sensors, spatial base schemes 
and nonlinear dissipation models are standalone modules. Therefore, a whole class of low 
dissipative high order schemes can be derived with ease. Unlike standard shock-capturing 
and/or hybrid shock-capturing methods, the nonlinear filter method requires one Riemann 
solve per dimension, independent of time discretizations. The nonlinear filter method is 
more efficient than its shock-capturing method counterparts employing the same order of 
the respective methods. 

Recently, these filter schemes were proven to be well-balanced schemes [59] in the 
sense that these schemes preserve exactly certain steady state solutions of the chemical 
nonequilibrium governing equation. With this added property these filter schemes can 
better minimize spurious numerics in reacting flows containing mixed steady shocks and 
unsteady turbulence with shocklet components than standard non-well-balanced shock- 
capturing schemes. In addition, for some stiff reacting flow test cases, the high order filter 
scheme is able to obtain the correct propagation speed of discontinuities whereas the stan- 
dard high order WENO scheme cannot [27, 74]. 

For simplicity of the presentation the discussion for the base scheme and post- 
processing step of the filter scheme is restricted to the inviscid part of the Navier-Stokes 
equations. For viscous gas dynamics the same order of spatial centered base scheme for the 
convection terms and the viscous terms are employed. For all of the LES computations the 
classical fourth-order Runge-Kutta time discretization is employed. 


Consider the 3-D compressible Euler equations in Cartesian geometry, 


U t + V • F = 0; U 



( pu \ 

puu T + p . 
\ u (e + p) J 


( 1 ) 


Here the velocity vector u = (u, v, w) T , the momentum vector m = (pu, pv, pw), p is the 
density, and e is the total energy. 

In a Cartesian grid denote the grid indices for the three spatial directions as (j, k, l ). 
The spatial base scheme to approximate the x inviscid flux derivatives F(U) X (with the grid 
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indices k and l for the y- and ^-directions suppressed) is written as, e.g.. 


OF 

dx 


DosFj, 


( 2 ) 


where D 0 8 is the standard eighth-order accurate centered difference operator. See [54] for 
the split form of 2. 

After the completion of a full Runge-Kutta time step of the base scheme step, the second 
step is to adaptively apply a post-processing nonlinear filter. The nonlinear filter can be 
obtained e.g., in the ^-direction by taking the full seventh-order WENO scheme (WEN07) 
[24] for the inviscid flux derivative in the £ -direction and subtracting D 08 Fj. The final 
update of the solution is (with the numerical fluxes in the y- and ^-directions suppressed as 
well as their corresponding y- and ^-directions indices on the x inviscid flux suppressed) 


At 


= u; M - —{h - +1!2 - 


(?) 


The nonlinear filter numerical fluxes usually involve the use of field-by-field approxi- 
mate Riemann solvers. If the Roe type of approximate Riemann solver [43] is employed, 
for example, the x-filter numerical flux vector H J+ \/ 2 evaluated at the U* solution from the 
base scheme step is 


Hj+ 1/2 = Rj+l/2Hj+l/2, 

where Rj+ 1/2 is the matrix of right eigenvectors of the Jacobian of the inviscid flux vector in 
terms of the Roe’s average states. Denote the elements of the vector H j+ 1 / 2 by //■ + , /2 , l = 
1, 2, ..., 5. The nonlinear portion of the filter li- + 1 / 2 has the form 

hj+l/2 = 2 UJ j+l/2 < fij+l/2- (4) 

Here u l j+1 , 2 is the wavelet flow sensor to activate the nonlinear numerical dissipation 
\<t> l j + 1/2 and the original formulation for n is a positive parameter that is less than or equal 
to one. Some tuning of the parameter k is needed for different flow types. A local n to be 
discussed next, depending on the local Mach number for low speed flows and depending 
on local shock strength for high speed flows, would minimize the tuning of parameters. A 
local flow sensor was discussed by Lo et al. [32] by taking advantage of the Ducros et al. 
shock flow sensor [1 1] to obtain a local artificial compression method (ACM) sensor for the 
original Yee et al. filter scheme [63]. 

The dissipative portion of the nonlinear filter \4>^ +1 , 2 = r/ / , /2 — // ^ , /2 is the dissipative 
portion of, e.g., WEN07 for the local /th-characteristic wave. Here fj l j + \/ 2 ar| d ^j+ 1/2 are 
numerical fluxes of WEN07 and the eighth-order central scheme for the Zth characteristic, 
respectively. Hereafter, we denote this filter scheme as WEN07fi. 
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Mach sensors, e=0 



Figure 3: Mach number sensors. f(M) (blue) function by Li and Gu, fi(M) (red) modified /(M), and 
/ 2 (M) (black) (includes low supersonic Mach numbers). 


A summary of the three basic steps to obtain o/ +]/ , 2 can be found in [51, 68]. For ex- 
ample, the flow sensor u;*- +1 , 2 to activate the shock-capturing dissipation using the cut off 
procedure is a vector (if applied dimension-by-dimension) consisting of “l’s” and “0’s”. 
For all of the computations, a three-level second-order Harten multiresolution wavelet de- 
composition of the computed density and pressure is used as the flow sensor [51]. 

2.2. Improved High Order Filter Method 

Previous numerical experiments on a wide range of flow conditions [63, 64, 51, 67, 68] 
indicated that the original filter scheme improves the overall accuracy of the computation 
compared with standard shock-capturing schemes of the same order. Studies found that 
the improved accuracy is more pronounced if the parameter k in (4) is tuned according 
to the flow type locally. For hypersonic flows with strong shocks, k is set to 1. For high 
subsonic and supersonic flows with strong shocks, k is in the range of (0.3, 0.9). For low 
speed turbulent flows without shocks or long time integration of smooth flows, k can be 
one to two orders of magnitude smaller than 1. In other words, k should be flow location, 
shock strength and local flow type dependent. The improved k proposed in [70] consists of 
a simple global k for smooth flows and a local k for problems with shocks and turbulence. 

2.2.1. An efficient global nfor low Mach number and smooth flows: 

The flow speed indicator formula of Li & Gu to overcome the shortcomings of “low 
speed Roe scheme” [30] was modified to obtain an improved global k denoted by k for (4) 
to minimize the tuning of the original n for low Mach number flows. 7c has the form: 


8 


with 


« = /i(M)k, 


(5) 


/i(M) 


min 


/M 2 y/4 + (1 - M 2 ) 2 
1 + M 2 


(6) 


Here M is the maximum Mach number of the entire computational domain at each stage 
of the time evolution. fi(M) has the same form as [30] except there is an extra factor “4p’ 
added to the first argument on the right-hand-side of the original form f(M) in equation 
(18) of [30]. The added factor provides a similar value of the tuning k observed from nu- 
merical experimentation reported in aforementioned cited references. With the flow speed 
indicator f\ (M) in front of k, the same n used for the supersonic shock problem can be used 
without any tuning for the very low speed turbulent flow cases. Another minor modification 
of the above is, 


fi(M) = max 


min 


/M 2 yfl + (1 — M 2 ) 2 
\~2 1 + M 2 


5 


where e is a small threshold value to avoid completely switching off the dissipation. A 
function which retains the majority of f\ (M) but includes larger Mach number for not very 
strong shocks is 

/ 2 (M) = (Q(M,2) + Q(M,3.5))/2 


or 

h W) = ma x((Q(M, 2) + Q(M, 3.5))/2, e), 


where 


The polynomial 


Q(M, a ) 


P(M/a) M < a 

1 otherwise 


P(x) = x 4 (35 - 84a; + 70a; 2 - 20a; 3 ) 


is monotonically increasing from P( 0) = 0 to P( 1) = 1 and has the property that P'(x) 
has three continuous derivatives at x = 0 and at x = 1. 


Below supersonic speeds, a simple and efficient global k can be obtained according to 
the maximum Mach number of the entire flow field and the value is determined by f\ (M) 
or / 2 (M) for non-zero a/ , |/2 . It is noted that if the original f(M) were used instead of 
fi(M) or / 2 (M) in Eq.(5), the amount of nonlinear filter dissipation could be too large for 
very low speed turbulent flows (for the same fixed k). See Fig. 3 for details. 
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2.2.2. Local flow sensor for a wide spectrum of flow speed and shock strength 

At each time step and grid point, the aforementioned global 7c is not sufficient to reduce 
the amount of numerical dissipation where needed for flows that contains a variety of flow 
features. A more appropriate approach is to obtain a“local k ” that is determined according 
to the above at each grid point. If known, a dominating shock jump variable should be used 
for shock detections. In other words, the filter numerical flux indicated in Eq.(4) is replaced 
by: 


h j+ 1/2 


9 K + 1/2^ + 1/2^+1/2]- 


(7) 


In the case of unknown physics and without experimental data or theory for compari- 
son, k l j + i/ 2 has to depend on the local Mach number in low speed or smooth flow regions, 
depend on local shock strength in shock regions and depend on turbulent fluctuations in 
vortical regions in order to minimize the tuning of parameters. According to the flow type 
locally, for each non-zero wavelet indicator uJ l j+1 / 2 , ^+ 1/2 should provide the aforemen- 
tioned amount (between (0, 1)) to be filtered by the shock-capturing dissipation <t> l j+1 / 2 - F° r 
problems containing turbulence and strong shocks, the shock strength should come into 
play. One measure of the shock strength can be based on the numerical Schlieren formula 
[21] for the chosen variables that exhibit the strongest shock strength. In the vicinity of 
turbulent fluctuation locations, n l - +l , 2 will be kept to the same order as in the nearly incom- 
pressible case, except in the vicinity of high shear and shocklets. 


Due to the fact that k works well for local Mach number below 0.4, k only needs to be 
modified in regions that are above 0.4. In other words, the final value of n l - +l , 2 is determined 
by the previous local 77 , if the local Mach number is below 0.4. If the local Mach number 
is above 0.4, at discontinuities detected by the non-zero wavelet indicator uj l j +1 , 2 , k 7 +1/2 
is determined by the shock strength (normalized between (0, 1)) based on the Schlieren 
formula near discontinuities. Again, if known, dominating shock jump variables should be 
used for shock detections. At locations with turbulence, determined by the turbulent sensor 
(e.g., u l ]+l j 2 obtained from employing wavelets with higher order vanishing moments), 
K l j + 1 / 2 is kept to the same order as in the nearly incompressible case, except in the vicinity 
of high shear and shocklet locations, where a slightly larger n l j+ i/ 2 would be used. Methods 
in detecting turbulent flow can be (a) Wavelets with higher order vanishing moments, and 
(b) Wavelet based Coherent Vortex Extraction (CVE) of Farge et al. [15] (Split the flow 
into two parts: Active coherent vortices and incoherent background flow). 

Results by the local k 1 j +1 , 2 that take the local flow speed and shock strength into consid- 
eration will be reported in [72], an expanded version of [70]. Preliminary study with more 
complex shock turbulence problems and the applicability of even wider flow types indicates 
the necessity of the local k l j +1 / 2 ■ 
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3. LES implementation 

In terms of turbulence modeling, there has been considerable progress in the develop- 
ment and usage of LES for the simulation of turbulent flows in the past few decades. This 
has been facilitated by the substantial increase in computing power. The triumphant jour- 
ney of LES started with the pioneering work of Smagorinsky [55], Lilly [31], Deardoff [10], 
Germano [19] and others. Comprehensive accounts on LES are provided by Sagaut [46] 
and Pope [35, 42] and reviews at different stages of the development are provided by Ro- 
gallo and Moin [44], Galperin and Orszag [18], Lesieur and Metais [29], and Meneveau and 
Katz [34]. The LES model used for the current simulation is discussed here with subgrid 
model summaries in Appendix A. 

In LES, the large-scale field can be obtained from the solution of the filtered Navier- 
Stokes equations, whereas scales smaller than the grid size are modeled. The filtering 
operation, which defines the large-scale variables (denoted by overbar), can be written as 





V, A) <p{y) dy , 


where D is the flow domain, G is the filter function, and A is the filter-width associated 
with the wavelength of the smallest scale retained by the filtering operation. 

Lor compressible flows, in order to avoid unclosed SGS terms in the continuity equation, 
a density- weighted Lavre filter operator is introduced. This operator is defined as tp = Jxp/ p. 

Lavre-filtered continuity, momentum and total energy equations in conservative form 
are respectively, 


dp_ dpuj 
dt dxi 


( 8 ) 


d (puj) d(pujUj) dp_ 

dt dxj dxi 


ddjj _ d Tjj 
dxj dxj 


(9) 


d(pE) d 


dt 


+ 


dx 


(, pE + p)uj - dijUi + q 3 


3 u 


d 


dx 3 




( 10 ) 


Unlike the ’bar’ and the ’tilde’, the ’breve’ -symbol does not denote a filter operation but 
indicates that the quantity is based on primitive filtered variables. Thus E refers to the 
resolved total energy, which is not equal to the filtered total energy. The resolved viscous 
stress tensor d VJ and the heat flux are defined as 


dij = 2 p(T) 


S, 


*? 



with Sij = (djUi + djUj) /2 


and 


dj 


-A (T) 


dT 

dxj 
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where p(T) and A (T) are the viscosity and thermal conductivity corresponding to the fil- 
tered temperature T. The constitutive equations are 


pE 



+ -pujUj and p 


Cv { 7 - 1 )pT. 


where heat capacity at constant volume, c v , and the ratio of heat capacities, 7 , are given 
constants. The subgrid scale stress is defined as 


Tij = p(UiUj - UiUj ), (11) 

and the subgrid terms of the energy equation are 

Qj = p^c v (ujT - UjT^J ( 12 ) 

Jj = Qj ~ Qj ( 13 ) 

P] = a jk u k - aJkUk (14) 

Sj = p {upu^k ~ UjUkUk) ■ ( 15 ) 


We will model the subgrid terms in order to close the system. The model for the subgrid 
stress, t^, employs the eddy viscosity hypothesis with an eddy viscosity, //,, with either a 
simple Smagorinsky model or a more advanced dynamic model. The details are given in 
Appendix A. In the energy equation, the subgrid terms J 3 , and £ 3 are set equal to zero, 
and the subgrid scale heat flux is modeled using the eddy diffusivity hypothesis, 

_ pgc v dT 
j Pr t dx 3 ’ 

where the turbulent Prandtl number, Pr t , is a given constant. 


4. LES of Temporally Evolving Compressible Mixing Layers (TML) 

In this section, employing WEN07fi, LES of temporally evolving mixing layers be- 
tween two streams moving with opposite velocities is considered, with U\ = —U 2 = A U /2. 
From here on, WEN07fi refers to the pre-processing step (Ducros et al. splitting of the in- 
viscid flux derivative) in conjunction with the eighth-order central spatial base scheme with 
the dissipative portion of WEN07 and the global flow sensor discussed in Section 2 with 
k — 0.7. The three main characteristics of compressible mixing layers are: 1) the self sim- 
ilarity property, which is characterized by linear growth of the layer width as well as the 
mean velocities and turbulent statistics being independent of the downstream distance when 
normalized by appropriate length and velocity scales, 2 ) the compressibility effects through 
the turbulence damping and the decrease of the mixing-layer growth rate for high convective 
Mach numbers, and 3) the presence of large-scale structure with shocklets. These organized 
structures play an important role in the dynamics of the mixing layer, its spreading and en- 
ergy transport. The objective of the current investigation is to verify that WEN07fi used in 
this study is capable of capturing the three-key points cited above. 
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4.1. Problem Setup 

The configuration of the TML is shown in Fig. 4. Five test cases (denoted LES-Q, 
i = 1, 5) are carried out with different convective Mach numbers (M c ) ranging from the 
incompressible case M c = 0.1 up to the supersonic one M c = 1.5. The latter corresponds 
to a highly compressible mixing layers, whereas the first two cases M c = 0.1 and M c = 0.3 
can be considered as quasi-incompressible and are used to compare with the experimental 
results of an incompressible shear layer. All of the simulations described below are per- 
formed at an initial Reynolds number, i?e Wo , based on the mean velocity difference A U, the 
average viscosity of the free streams and the vorticity thickness 5 Uo of 800 with 8 Wo =4 5q 0 , 
where = AU / ( du/dy) max is the vorticity thickness of the shear layer, and do is the mo- 
mentum thickness given by Eq. 17 (see later section). It is worth noticing that Re u reaches 
values as large as 3 x 10 5 at the end of the simulation, which is one order of magnitude 
higher than similar DNS and LES computations reported in the literature [39, 33, 16]. Ta- 
ble 1 summarizes the detail of flow parameters for both LES and previous DNS data of the 
literature. The mean flow is initialized with a tangent hyperbolic profile for the streamwise 
velocity, u(y) = \ AU tanh [y/( 2 Aj], while the two other velocity components are set to 
zero. In addition to these mean values, three-dimensional turbulent fluctuations (i//, v' . w') 
are imposed, while initial pressure and density are set constant. Since the simulation is tem- 
poral, the initial velocity perturbations are computed using a digital filter technique [26]. 
This procedure utilizes the prescribed Reynolds stress tensor and length scales of the prob- 
lem concerned to generate the corresponding fluctuating velocity field, taking into account 
the nature of the autocorrelation function for the prevailing turbulence. The digital filter 
algorithm is given in Appendix B. The length scales are chosen as 5 W0 in each direction. 
The Reynolds stresses have a Gaussian shape in y with amplitudes chosen to be similar to 
the experimental peak intensities observed in incompressible mixing layer [5]. 

Periodic boundary conditions are enforced in the streamwise (x) and spanwise (z) di- 
rections, while non-reflecting conditions are applied at both top and bottom boundaries (y 
direction). The use of a periodic boundary condition in the x direction corresponds to the 
temporal formulation of mixing layer evolution, which evolves only in time as it spreads in 
V- 

4.2. Mesh Requirements 

Similar to [16], a computational domain of lengths L x x L y x L z = 1200 5g 0 x 370 S# 0 x 
270 So 0 is used with the corresponding mesh points N x x N y x N z — 512 x 211 x 131. 
The same grid system uniformly spaced in the x and z directions and stretched in the y 
direction is employed for all considered cases. The stretching function for the ^/-direction 
is based on where L y is the box size in the ^/-direction and the stretching factor 

b y = 3.4. The mapped coordinate q is equally spaced and runs from —1 to 1. The high 
resolution (HR) grid used in this study contains an order of magnitude fewer cells than that 
of the DNS of Pantano and Sarkar [39] compared to the domain length. The emphasis of 
the HR simulation is to produce an LES solution that predicts the trends of the DNS as well 
as experimental data for both quasi-incompressible and highly compressible mixing layers. 
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Figure 4: Schematic configuration of temporal mixing layer 

To ensure that the computational domain in the x- and ^-directions is sufficiently wide, the 
two-point correlation functions are analyzed 

N-k r 

Rv>cp(r) = VkV'k+k^ Ay = 0,1,...., TV- 1, (16) 

k= 1 

where r = k r A, N is half of the number of grid points in the homogeneous directions with 
grid size A and <p' represents the fluctuations of flow variables. 

The computed two-point autocorrelation coefficients R w {r)/R w ( 0) (pressure as well 
as velocity components) in the homogeneous directions (x and z) are reported in Fig. 5 
as a function of the distance in the stream- and spanwise coordinates at the middle of the 
mixing layer (L y / 2). This figure shows that the flow variables are sufficiently decorrelated 
over distances L x /2 and L z / 2, thus ensuring that the streamwise as well as the spanwise 
extents of the computational domain are sufficient to not inhibit turbulence dynamics. Also, 
the length in the y direction is selected to be large enough for the flow to achieve a fully 
developed state before the effects of the upper and lower boundaries are felt. In terms of 
turbulent length scales, the Kolmogorov length scale r) and an average (isotropic) Taylor 
micro-scale A are defined as y = (z/ 3 /^) 1 ^ 1 and A = yA5 v k/e, where k = \{u' 2 + v' 2 + 
w' 2 ) is the turbulent kinetic energy. The computed integral length scales (A.,. . A z ) and the 
Kolmogorov scale are also given in Table 1 for further comparison. In our case the integral 
scales are given by 

A _ [ Lx/2 (Ui(Xi,t)Ui(Xi +p,t)) A A _ f Lz/2 (Ui(Xi,t) Ui(Xi +P,t)} J 

./o M) dP ’ Z ~ Jo (Uf) dp ’ 

where p is the distance between two points in the flow. The integral length scale is impor- 
tant in characterizing the structure of turbulence. It is a measure of the longest correlation 


14 


Table 1 : Flow parameters and turbulent length scales during the quasi-self-similar stage. 


Case AI C R&u 0 L x / A. x L z /A. z Aij tnxn f A Ay mzn / t\ 


DNS [45] 

0 

800 

- 

DNS [39] 

0.3 

640 

20.4 

LES -Cl 

0.1 

800 

30.2 

LES-C2 

0.3 

800 

30.2 

LES-C3 

0.8 

800 

12.1 

LES-C4 

1.0 

800 

12.0 

LES-C5 

1.5 

800 

12.0 


15.3 

> 1 

« 1 

15.1 

0.64 

49.3 

15.1 

0.65 

49.5 

11.0 

0.74 

50.1 

11.4 

0.67 

53.7 

10.2 

0.82 

55.2 


distance between the flow velocity (or vorticity, etc) at two points in the flow field. Recent 
work concludes that a reasonable lower limit on the domain must be at least six times larger 
than the integral length [37]. This recommendation is consistent with the data shown in 
Table 1, where the spatial domain is between ten and thirty times larger than the integral 
length. Also, it is evident from Table 1 that the integral lengths are sufficiently small com- 
pared to the computational domain and the grid resolution is adequate to resolve the large 
scales of turbulence. 

Owing to the high computational cost of the simulations, the numerical code is fully 
parallelized running on up to 600 processors. In total, the present simulation required 2000 
CPU hours each on modem SGI IC, Pleiades and Columbia supercomputers at NASA Ames 
Research Center. 

4.3. Mean Flow and Turbulent Statistics 

LES computations are carried out up to dimensionless time r = tAU /Sg 0 ~ 3000 for 
the higher convective Mach number cases and r ~ 1200 for the quasi-incompressible cases. 
In order to compare the LES results with experimental data, the time averaged flow quan- 
tities (</?) and (Ip) are extracted from the flow field during the self- similar time period 
(600 < r < 1000 for LES-C1 and LES-C2 and 2000 < r < 2800 for LES-C3, LES- 
C4 and LES-C5). Note that throughout this paper only resolved quantities are considered; 
subgrid-scale contributions are not added onto e.g. the turbulent stresses. To validate the 
low-Mach-number LES case, previous DNS studies of the incompressible shear layer, in- 
cluding Rogers and Moser [45], Pantano and Sarkar [39], as well as experimental studies 
by Bell and Mehta [5] and Spencer and Jones [56], are used. Further experimental results 
on the compressible shear layer (Papamoschou and Roshko [40], Elliot and Samimy [13], 
Barre et al. [4], Chambres et al. [8]) and DNS results obtained by Pantano and Sarkar [39] 
are used to compare with the high-Mach-number simulations. 

As recommended by Rogers and Moser [45], the momentum thickness Sg is used for 
self-similar scaling rather than the vorticity thickness 8 U , because it is less sensitive to sta- 
tistical noise as it is an integral quantity evolving smoothly in time, while 8 U is obtained 
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Figure 6: Time evolution of normalized momentum thickness for LES-C2 compared to the DNS of Pantano 
& Sarkar [39] for M c = 0.3. 

from the derivative of the mean velocity and may exhibit oscillations during the flow evolu- 
tion. Therefore, the time evolution of the momentum thickness of the flow calculated using 
the definition 


6e = 





(17) 


is shown in Fig. 6. Excellent agreement with the DNS simulation is obtained and the linear 
slope is recovered after a short transient, showing the self-similar state of the mixing layer. 
The growth rate d(S 9 / S do ) / dr = 0.016 (slope of the linear curve fit). The ratio of vorticity 
thickness to momentum thickness is = 5^ /So — 4.5 with Re u = 603391 at r maa; ~ 
1200. This is in excellent agreement with the DNS growth rate of quasi-incompressible 
case with M c = 0.3 of Pantano & Sarkar [39]. Note that in Eq. 17, ( ) re j represents the 
reference state which is the arithmetic mean of the free stream 1 and 2. 

The normalized mean streamwise velocity for LES-C1 and LES-C2 are presented in 
Fig. 7. It can be seen that the two LES profiles collapse together and are in excellent 
agreement with both DNS [39] and experimental results [5, 56]. 
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Further validation of the present LES is achieved through the comparison of turbulent 
intensities in the self-similar region (calculated by averaging over profiles plotted in simi- 
larity coordinates) seen in Fig. 8, where different components of the normalized Reynolds 
stress tensor, ^ Iff/ All ( R, :J = pu"u" /p), are compared to DNS and experimental data. 
Interestingly, at Me = 0.1 the LES agrees better with experimental [5, 56] and DNS [45] 
data in the incompressible shear layer than the DNS results of Pantano & Sarkar [39]. Ta- 
ble 2 summarizes the comparison of peak turbulent intensities, as well as the anisotropic 
deviation on the centerline of the layer. It is evident that very good agreement between the 
present LES and previous results is obtained for this measure of anisotropy. It is important 
to notice that both LES-C1 and LES-C2 give almost the same results, probably because 
both are in the incompressible (or weakly-compressible) regime. 



Figure 7: Comparison of normalized mean streamwise velocity for LES-C1 and LES-C2. 


4.4. Compressibility Effects 

Apart from studying the self- similarity property in the mixing layer, effects of con- 
vective Mach number are also investigated using LES. Fig. 9 shows the time evolution 
of momentum thickness for various convective Mach numbers. It can be seen that af- 
ter a relatively long time (r > 2000 compared to the incompressible one), correspond- 
ing to the initial transient, the mixing layer grows quasi-linearly with spread rates of 
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a/OWA U y/{R^}/AU 


Table 2: Comparison of peak turbulent intensities of incompressible mixing layer 



Bell et al. 

Pantano et al. 

LES-C1 / LES-C2 

M c 

0 

0.3 

0.1/ 0.3 

^(R^/AU 

0.18 

0.155 

0.17 

yf{R^)/AU 

0.14 

0.134 

0.134 

V(R^)/AU 

0.146 

0.143 

0.143 

^/{Ruj/AU 

0.10 

0.103 

0.106 

yj (-R 22 )/ (-R 11 ) 

0.777 

0.788 

0.788 


0.555 

0.606 

0.623 
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Figure 8: Comparison of the normalized Reynolds stress tensor, = (pu'u'} / (p) . 
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d{S e /5e 0 )/dr = 0 . 0165 , 0 . 0101 , 0.0084 and 0.0075 for cases LES-C2, FES-C3, LES-C4 
and LES-C5, respectively. 

In compressible mixing layers, all of the assessments of compressibility effects can 
be related to the convective Mach number, M c , through the compressibility factor, <f> = 
(dS /dr) c / (dS /dr)i, which is the ratio of the compressible growth rate to the incompress- 
ible growth rate at the same velocity and temperature ratios. The calculated compressibility 
factor is significantly less than one (the incompressible counterpart), the ratio of the two 
being less than 0.43 for M c > 1 . This is consistent with previous findings on the effects of 
compressibility on mixing-layer growth rate such as the nonlinear regression fit of Barone 
[3, 28] plotted in Fig. 10. This plot shows the ratio of compressible mixing-layer growth 
to incompressible mixing-layer growth rate as a function of M c , and data from different 
experiments and previous DNS have been included for comparison. The three higher com- 
pressibility cases have growth rates that agree well with previously published data. As 
already pointed out by Papamoschou [41], the growth-rate reduction starts at subsonic val- 
ues of M c and is evidently completed before M c becomes supersonic (Fig. 10). This implies 
that compressibility takes effect before any shock or expansion waves appear in the flow, in 
the convective frame of reference. 

It is worth noticing that the data in this figure exhibit significant scattering that is partly 
attributable to the different experimental conditions. As pointed out by Barone et al. [3], 
one can mention that future investigations should be conducted at higher convective Mach 
numbers to better determine the asymptotic value of ( I>. 

Also, the dependence of the turbulent kinetic energy of the shear layer on M c is shown 
in Fig. 11. The simulations show that the turbulence intensity decreases with increasing 
convective Mach number. The decreased level of energy is responsible for the reduction 
of the mixing thickness growth rate as already pointed out by Samimy [49], Vreman et al. 
[58] and many other FES and DNS studies [39, 16, 50]. 


4.5. Flow Structures and Shocklets 

The invariant of velocity gradient tensor Q and the corresponding normalized form A 
are defined by 


Q 2 I '3) j S ; ] , A 




[QijQij + SijSij] 

where S{j T Aa ) ^0,0 /^* 


(18) 


The iso-surfaces of Q and A are plotted for flow visualization of mixing layers. It is evident 
that the positive values of Q and A represent the vortex dominated portion of the flow. 
Three-dimensional perspective views of iso-surfaces of Q are presented in Fig. 12 for FES- 
C2 in a self-similar state. The 3D complex vortex tube structures are clearly evident from 
these figures. 

With regard to the highly compressible case M c = 1.5, the complexity of three- 
dimensional flow structure leads to difficulties in the identification of shocklets. One good 
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Figure 9: Time evolution of normalized momentum thickness for LES at different convective numbers. 
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Compressibility factor 



Figure 10: Temporal mixing: Compressibility factor as function of the convective Mach number from different 
experimental mixing-layer studies selected by Barone et al. [3]: (a) Bogdanoff [7]; Papamoschou and Roshko 
[40]; (b) Chinzei at al. [9]; (c) Samimy and Elliott [48, 49]; — nonlinear regression curve from [3] with 
4>(M C ) = 1 - <zi [1 - 1/(1 + a 2 M“3)], ai = 0.5537, a 2 = 31.79, a 3 = 8.426. fd) Gruber et al. [20], LES 
computations by WEN07fi (red solid circles) for M c = 0.1, 0.3, 0.8, 1.0, 1.5. 
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Figure 11: Normalized turbulent kinetic energy at different convective Mach numbers. 
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method to identify the location of a shock is to use a Schlieren based technique to portray 
shocks and even more weak discontinuities in the fluid (see ref. [21] for more details about 
flow visualization). Since, in our case, the initial density is uniform, we used Schlieren 
based on the dilatation of the velocity field Vu to highlight the eddy shocklets (see Fig. 
14). Note that shocklets start to appear at Mach number less than unity, i.e. in the lower 
part of the transonic regime. As expected, for higher Mach numbers the shocklets become 
stronger and are preferentially organized in oblique waves (see Fig. 14), corresponding to 
stationary inviscid shocks at a dominant propagation directions, 9, and a nominal Mach 
number, M n = A U / (2c), where c is the speed of sound in the unperturbed region. Avital 
et al. [1, 2] provide a correlation to find the propagation angle corresponding to the most 
perturbed waves. It is given by M c cos9 = 0.6. Therefore for M c = 1.5 — > 9 t h = 66°. In 
our case, a visual measurement of the oblique wave angle gives an approximate value of 
9 sim ~ 65°, which is very close to the predicted value. 

From the present computation it can be seen that oblique structures start to occur at con- 
vective Mach numbers less than unity. These structures are related to compression waves 
emanating from the shear layer and also the existence of other perturbing pressure distur- 
bances and lead to enhanced mixing through the creation of streamwise vortices. 

5. Concluding Remarks 

This paper illustrates some recent progress in computations of compressible turbulence 
using a high-order spatial scheme on a LES model for temporally evolving turbulent mixing 
layers. Results obtained including flow visualization, streamwise velocities, fluctuating 
velocities and Reynolds stresses agree well with experimental results. The current LES 
agree with the previous DNS of the mixing layer by Vreman et al. [58] and Freund et al. 
[17] that show decreased turbulence production with increasing M c . 

The present study serves as a validation and performance of the improved filter schemes 
of [70] on a representative complex compressible turbulent flow consisting of a wide range 
of flow speeds. All the computations use the Ducros et al. splitting of the inviscid flux 
derivatives and WEN07fi with k and n = 0.7 described in Section 2.2.1. In all M c cases, no 
tuning of WEN07A scheme parameters were needed. LES comparison among WEN07A, 
WEN05 and WEN07 for the TML is reported in [71]. Studies indicated that WEN07A 
compared well with experimental data and published DNS work. For all the considered 
M c cases, solutions by WEN05 and WEN07 compared poorly with experimental data and 
DNS computations. The comparative study among WEN07fi, WEN05 and WEN07 is 
reported in [71] was the first step in determining the suitable order of filter schemes to be 
used for the current physics based study. 

The same high order filter scheme is being used for the simulation of two much higher 
M c cases of M c = 2, 3. The computational box size, especially in the y-dircction has to 
be doubled or more. A finer grid is also needed in order to obtain an accurate and stable 
solution. These computations are many times more CPU-intensive than the lower M c cases. 
Results will be reported in a forthcoming paper. 
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Side view 



Top view 



Lateral view 



Figure 12: Iso-surfaces of Q = 0.01, Q max at r = 1000, LES-C2. 
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Figure 13: Instantaneous dilatation flow-field at r = 1000 for three different convective Mach numbers 
(M c = 1.5, 1.0, and 0.8 from top to bottom). Note that the plot are based on the non-linear dimensionless 
variable, </> = 1 — tanh [k V- u/(V- u) ma J, the parameter k = 0.5 governs the amplification of small 
gradients; a value of k close to 15 provides good results. 


26 



wave angle 


Eddy shocklets 


Turbulent mixing 


Figure 14: Instantaneous numerical Schlieren pictures at r = 2000, LES-C5. 
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Appendix A. Subgrid Model 


The most widely used and simplest model is Smagorinsky’s model [55], which employs 
an eddy viscosity hypothesis to express the subgrid scale stress as, 


Tsj r^kk^ij 2 /// ( Sij ^ S kk S tJ ) . . 

The eddy viscosity, p t , is modeled according to, 

Pt = pC s A 2 \S\, 

where C s is the Smagorinsky constant and \S\ is defined as: 

1/2 


(19) 


\S\ = 2 SijS, 


The model for the isotropic part of the subgrid scale stress was proposed by Yoshizawa 
[75] as 


r kk = 2C I pA 2 \S\ 2 . 


According to Erlebacher et al. [14], r^k can be neglected in flows where turbulent Mach 
number, M t = is less than 0.4. 

The dynamic procedures have been developed to evaluate the parameters used in the 
subgrid models using the information provided by the resolved scales. The original proce- 
dure was developed by Germano et al. [19] and later modified by Lilly [3 1 ] for incompress- 
ible flows. Moin et al. [36] generalized the dynamic procedure for compressible flows. The 
procedure is as follows: The sub grid scale stress for compressible flows which is defined in 
Eq. 1 1 can be rewritten as 


T{j pUiUj 


ptl, pUi 

p 


( 20 ) 


Now the field Ui is considered as instantaneous field and a test filter with the filter kernel 
G(x — y, 2A) is applied to the LES filtered Navier-Stokes equation to get the resolved 
turbulent stress, 





pili pUj 

p 


( 21 ) 
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where q denotes the test filtered variable q. The subtest stresses are defined by, 

Tij pu.iUj 


pUi pUj 


p 

From Eqs. 20, 22 and 21, we get the Germano identity, 


T —T 

— j-ij 


T, 




( 22 ) 


(23) 


In the above equation, the two terms on the r.h.s. can be modeled according to 
Smagorinsky. The term in the l.h.s. can be explicitly calculated by applying the test fil- 
ter to the simulation results which were obtained using the first filter. 

The anisotropic part of the subgrid stress given in Eq. 22 can be modeled according to 
Smagorinsky as 

Tij - -TiAj = —2C s J)A 2 \S\S*j (24) 


with S'* = Sij — \Su8ij. The isotropic part of the subgrid stress given in Eq. 22 can be 
modeled according to Yoshisawa as 


T u = 2C/pA 2 |S| L . 


(25) 


After applying the test filter to the subgrid stress tensor of tq (Eq. 19) and substituting 
it along with Eq. 24 into Eq. 23, we get, 


L ij L,,d,i 6 


ll^ij 


T Cs 


-2pA 2 |S|S* + 2A 2 ( I p|S|S y ) - ; (mSa ) 4 


V 




(26) 


M c . s 

IJ 


La — Cj 


2pA i |S'p — 2A J p\S\ 


(27) 


M, 


Cl 


The above equations can be rewritten in compact form as 


= C,ME% L u = CjM^' . 


v 


*? 


(28) 


Modifying the original approach by Germano to find C s and 67/ from the above equa- 
tions, Lilly [31] introduced least square method which gives 


67, = 


(^ M S s )h 


C r = 


( L ll ) 


H 


« 7 >rr' 


(29) 


In the above equations, as we can see, an averaging is done in the homogenous direction. 
This is to avoid excessively large local values of 67 s and 67/ which may destabilize the 
numerical simulation. 
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Appendix B. Algorithm for Turbulence Initialization [26] 

(a) Choose length scales in each direction L x = n x Ax, L y = n y Ay and L z = n z Az 

(b) Choose a filter width Nf a > 2 n a , a = x,y,z 

(c) Initialize and store random fields 1Z n with zero mean and unity variance of dimension 
[-N fx : N fx ,—N fy + l:N fy + N y , -N fz + 1 :N fz + N z ] , where N y x N z are the 
number of the mesh points in yz plane 

(d) Calculate filter coefficients b(i,j, k), where b UJ j ; = b t ■ bj ■ b k . 


Y nj 

A,j=-Nf U j 


1/2 


and b k = exp 


nk 2 
'2 n 2 


(e) Apply the filter operation for j = 1, , N y . k = 1, , N z 

Nf x Nf y N fz 

Ua(j,k)= 5Z b ( i 'J', k ') n *( i ',j+j', k + k') 

i’=-N fx j’=—N fy k’=-N fz 

which results in the two-dimensional arrays of spatially correlated data U. a 

(f) Perform the following coordinate transformation to get U' a (j,k) = « tJ U n (j, k), 
with prescribed Reynolds stress tensor 

/ (i?n) 1/2 0 0 \ 

a ij — I Rll/ttll (R22 — °2l) 1//2 0 I 

\ Rzi/^U (R32 — (l 2 l(l 3 l) / a 22 (R33 ~ O31 — a 32) 1//2 J 

(g) Calculate u a (j, k ) = u a (j, k ) + U' a ( j , k) for first (j, k) plane 

(h) Discard the first y, z plane of TZ a and shift the whole data 

n a (■ i,j,k ) := TZ a (i + l,j, k) 

(i) Generate new random numbers to fill the plane 7 Z a ( Nf x ,j , k) 

(j) Repeat the steps (e) to ( i ) for each mesh point in the x direction 
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